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ABSTRACT 

This paper is the third in a series presenting the results of direct numerical integrations of the 
Fokker-Planck equation for stars orbiting a supermassive black hole (SBH) at the center of a galaxy. 

The algorithm of Paper II included diffusion coefficients that described the effects of random (“classi¬ 
cal”) and correlated (“resonant”) relaxation. In this paper, the diffusion coefficients of Paper II have 
been generalized to account for the effects of “anomalous relaxation,” the qualitatively different way 
in which eccentric orbits evolve in the regime of rapid relativistic precession. Two functional forms 
for the anomalous diffusion coefficients are investigated, based on power-law or exponential modifica¬ 
tions of the resonant diffusion coefficients. The parameters defining the modified coefficients are first 
constrained by comparing the results of Fokker-Planck integrations with previously-published V-body 
integrations. Steady-state solutions are then obtained via the Fokker-Planck equation for models with 
properties similar to those of the Milky Way nucleus. Inclusion of anomalous relaxation leads to the 
formation of less prominent cores than in the case of resonant relaxation alone, due to the lengthening 
of diffusion timescales for eccentric orbits. Steady-state capture rates of stars by the SBH are found 
to always be less, by as much as an order of magnitude, than capture rates in the presence of resonant 
relaxation alone. 


1. INTRODUCTION 

Paper I of this series (|Merrittll20I5all described a numerical algorithm for integrating the Fokker-Planck equation 
for f{E,L,t), the phase-space density of stars orbiting a supermassive black hole (SBH) at the cente r of a galaxy; E 
and L are respectively the orbital energy and angular momentum per unit mass of a star. Paper H (|Merrittl 120151)1 1 
presented steady-state and time-dependent solutions for / based on diffusion coefficients that describe the effects of 
both “classical” (random) and “resonant” (correlated) relaxation; the latter becomes progressively more important 
relative to the former as one moves inside the gravitational influence sphere of the SBH. A single value for the stellar 
mass, m*, was assumed in both papers. 

This paper extends the tre atment of Paper H t o include the qualitatively different sort of evolution experienced by 
eccentric orbits near a SBH (jMerritt et al.ll201l[l . Such orbits undergo rapid apsidal precession due to the effects of 
general relativity (GR), at an orbit-averaged rate 


^ _ 3(GM.)^/^ 

dt (1 — e^) (? 


( 1 ) 


(“Schwarzschild precession”: iMerrittI 120131 Equation 4.205). Here a and e are the orbital semimajor axis and eccen¬ 
tricity respectively and w is the argument of periapsis. Precession described by Equation m affects the collective 
evolutio n in two, distinct ways, ( i) The “coherence time” is defined as the mean precession time for orbits at a given 
energy (jRauch fc Tremainelll996| l. Near the SBH coherence times become progressively shorter due to Schwarzschild 
precession. This effect was correctly accounted for in Papers I and H and in many earlier treatments of resonant 
relaxation, (ii) At any energy, sufficiently eccentric orbits undergo apsidal precession on a shorter timescale than other 
orbits of the same energy, due to the strong eccentricity dependence of Equation O- Such high-eccentricity orbits 
might be expected to evolve in a manner qualitatively different than described by the equations of resonant relaxation, 
since their orientation with respect to the torquing potential changes in a time short compared with the coherence 
time. 

Since it is the high-eccentricity orbits that are most amenable to capture by the SBH, Schwarzschild precession was 
recognized early on a s a potentially important mediating factor with re gard to rates of capture from tightly-bound 
orbits around a SBH (iHopman fc Alexandedl20d^ iMadigan et ani2011ll . __ 

The first, fully self-consistent investigation of orbital evolution in this regime (jMerritt et al.l[20Til l revealed a new phe¬ 
nomenon. Stars near the SBH undergo random walks in angular momentum due to resonant relaxation, but when their 
eccentricities reach a certain maximum value (depending on a), their trajectories “bounce,” returning after roughly 
one coherence time to lower values of e, where they continue to evolve under the influence of resonant relaxation. The 
locus of reflection in the (a, e) plane was termed the “Schwarzschild barrier” (SB) and an approximate analytic expres¬ 
sion for its location, L = Lsb(A), was derived. Subseq uent studies have confirmed this phenomenon using differen t 
integration schemes for the A^-body equations of motion (jBrem et al.ll20l4 iHamers, Portegies Zwart fc Merrft3l20I4fl . 
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A characteristic of motion near and below the SB {L < Lsb{E)) is that the apsidal precession time is short compared 
with the coherence time, and with the time over which resonant relaxation would be able to change L in the absence of 
the rapid precession. Angular momentum evolution in the region below the SB was called “anomalous relaxation” by 
iHamers. Portegies Zwart &: Merritt ()2014l l. This name reflects the fact the the evolution in this regime is qualitatively 
different than the evolution described by the equations of either classical or resonant relaxation. For instance: diffusion 
rates in this regime drop rapidly wi th decreasing L, and there is a net drift in the direction of increasing L. 

The direct iV-body integrations of iMerritt et al.l (1201111 showed that the SB is not completely impermeable, although 
captures by the SBH were found to occur at a rate that was about an order of magnitude lower than in simulations 
that omitted the first post-Newtonian (IPN) terms from the equations of motion; that is, the terms that generate 
Schwarzschild precession. Extrapolating the capture rates in those small-simulations to real galaxies is not straight¬ 
forward. One reason is the absence, in the A^-body models, of stars initially distant from the SBH that would diffuse 
inward and replace those lost to the SBH, thus establishing a steady state. Another reason was pointed out by Hamers 
et al. (2014). At sufficiently low L, anomalous diffusion rates can become so low that classical relaxation once again 
sets the timescale for angular momentum evolution. Using an approximate test-particle algorithm, Hamers et al. were 
able to simulate systems of much larger N and to cleanly delineate three regimes of angular momentum evolution, at 
energies for which the SB exists: 

1. Lsb{E) ^ L < 1 (resonant relaxation) 

2. Lnr(A) < L < Lsb{E) (anomalous relaxation) 

3. Lic{E) < L < Lnf.(A) (classical relaxation) 


(see their Figure 1). Here Tnr is the angular momentum at which Schwarzschild precession is so rapid that the torques 
driving resonant relaxation are almost completely ineffective at changing L, so th at classical relaxation dominates the 
evolut ion once more. Lie is the angular momentum at the edge of the loss cone. IHamers. Portegies Zwart fc Merritt 
(|2014ll derived an approximate expression for L^b.{E) and showed that in the simulations of iM^ritt et al.l ( 201lir 
classical relaxation dominated the evolution in L over much of the (a, e) plane, including even some regions with 
L > Lsb{E). They argued that this fact would complicate the extrapolation of the A^-body results to real galaxies. 

This paper, the third in a series, addresses these issues by incorporating into the Fokker-Planck algorithm expressions 
for the diffusion coefficients that account for anomalous relaxation. By integrating /(A, L) forward in time using these 
new diffusion coefficients, steady-state solutions are constructed that are valid fully into the “Schwarzschild” regime 
defined in Paper H - roughly an order of magnitude nearer to the SBH than the solutions of Paper H, or indeed any 
other published simulation. 

Section [2] reviews the numerical algorithm used here; further details are given in Papers I and H. Section [3] presents 
the functional forms adopted for the anomalous diffusion coefficients. Since there does not yet exist a good theory 
for orbital evolution in this regime, different parametrized forms for the diffusion coefficients are considered and 
constrained by comparison with previously-published simulations. Section|3]presents steady-state solutions for f{E, L) 
with parameters chosen to describe the nuclear cluster of the Milky Way; the results are compared with those of Paper 
H that did not incorporate anomalous relaxation. Section [5] discusses some implications of the results obtained here 
and ^sums up. 


2. METHOD 

As in Papers I and H, stars are assumed to have a single mass, m*, and to be close enough to the black hole (SBH) 
that the gravitational potential defining their unperturbed orbits is 


(i)(r) =- - = —'tp{r) (2) 

r 

with M, the SBH mass, assumed constant in time. Unperturbed orbits respect the two isolating integr als A, the 
energy per unit mass, and L, the angular momentum per unit mass. Following iCohn fc KulsrudI (|1978[1 these are 
replaced by £ and TZ where 


f =-A =-y-hV'C?’), ^=^5 


( 3 ) 


Lc{£) is the angular momentum of a circular orbit of energy £ so that 0 < 72. < 1. £ and 72 are related to the 
semimajor axis a and eccentricity e of the Kepler orbit via 


a = 


GM, 

2 £ 


= 1 - 72. 


( 4 ) 


Spin of the SBH is ignored. 
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The time dependence of the phase-space number density of stars, f{£,TZ), is described by the orbit-averaged Fokker- 
Planck equation 


d f d d 

—4‘s = Dee-^ + Dstz-^ + Dsf, —4>n = + D-nn^^ + -Dtc/ 


(5) 


with flux coefficients 


D, = -(AQ - + 1T((A«^) + ij^(A«AK) , 

Dn = -{An) - ^{ASAU) + ^±{A£An) + ^^{(Anf) , 

Dee = \{{A£f) ,Den = Dne = \{A£An) ,Dnn = \{{Anf) (6) 


and J = y/2'K^G^M,^£~^/'^ (|MerrittlI2013L 5.5.1). Quantities in ( ) are orbit-averaged diffusion coefficients. The 
functional forms of the diffusion coefficients are discussed below. 

Loss of stars into the SBH is controlled by the choice of ric, the radius of the physical loss sphere around the SBH, 
and by the conditions imposed on / at the loss-cone boundary, TZ = TZ\c(£)i defined as 


7^lc(f) =2^ ( 1 - ) , £<£n, £n = 


£\c 


2£ic 


GM, 

2nc 


(7) 


TZ\c is the normalized angular momentum of an orbit with (Newtonian) periapsis at rn- The 7?.-directed flux of stars 
across the loss-cone boundary is 

F{£) d£ = -J{£)(j)n{nx,) d£ = -J{£) <\)nA£) d£. (8) 

Two quantities that play important roles in angular momentum diffusion near the loss-cone boundary are 2?, 


V{£)^ 


{{An)\ 


27^ 


Dnn{£, 72.1c) 


T^—TZic 


TZic 


and qic, 


„ . = P{£M£) 

^ TZn{£) 


(9) 


( 10 ) 


is effectively an orbit-averaged, angular momentum relaxation time at energy £. The quantity qic measures the 
change in angular momentum per orbital period, compared with the size of the loss cone. The loss-cone boundary 
conditions adopted in all the integrations presented here were the “Cohn-Kulsrud boundary conditions” defined in 
Paper I. No attempt is made to solve for / inside the loss cone, i.e. at 72 < 72ic, since / does not satisfy Jeans’s 
theorem in this region. 

Solutions are obtained on a {N^ x N^) grid in (X, Z), where 


X = In 7? = In 


' L Y 

L,{£) 


Z = \n{l + /3£*)=ln{l + l3£/c^) . 


Integrations presented here used = 64 grid points. The code adopts units such that 


G = M. = c = 1 


( 11 ) 

( 12 ) 


allowing the results to be scaled to different masses of the SBH. Dimensionless parameters that must be specified 
before the start of an integration include In A and 0ic = ric/vg. 

In Paper II, the diffusion coefficients had the forms 


(Af) = (A£:)ck, { iA £ f ) = { iA £ f ) cK , { A £ An ) = (A£:A72)ck, 

(A72) = (A72)ck + (A72)rr, {{ Altf ) = ((A72)^)ck + ((A72)^)rr . (13) 

The subscript CK indicates that the diffusion coefficient is compute d as in iCohn fc KulsrudI (1197811 : their derivation 
was based on standard as sumptions about randomn ess of encounters (IRosenbluth et al.lll957ll . The subscript RR refers 
to “resonant relaxation” ([Ranch fc Tremainelll996|l . The resonant diffusion coefficients were expressed as 


(A72)rr = 2A{£) (1 - 272), ((A72)")rr = 4A(£)72 (1 - 72). 


(14) 
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The term containing the E dependence is 


A{a) = al 


Mi, ^ 1 tcoh 

W, 


as = 1 . 6 , 


a = 


GM, 


(15) 


Here N = N{r < a) is the number of stars instantaneously at radii than a, M* = 'rni,N, P is the Kepler (radial) 
period, and tcoh is the coherence time, defined as 


+-!=+-! _L+-i 

''coh — ‘•coh,M ^ ‘•coh,S 


M, 

^coh,M(^) — “77 ^ 

Nm^ 


^coh,s(^) — 


12 r„ 


-P. 


(16) 


tcoh,M is the mean precession time for stars of semimajor axis a due to the distributed mass around the SBH (“mass 
precession”), and tcoh,s is the mean precession time due to the IPN corrections to the Newtonian equations of motion 
(“Schwarzschild precession”). 


3. ANOMALOUS DIFFUSION COEFFICIENTS 

The diffusion coefficients m are affected by general relativity (GR) to the extent that GR determines the coherence 
time via equation (IMl). Another GR-related phenomenon is the Schwarzschild barrier (SB ), the tendency of orbits 
near the SBH to avoid high eccentricities. The SB was first observed in A^-body simulations (jMerritt et al. 1 1201111 . as a 
locus in the (E, L) plane where trajectories “bounced” during the course of their random walks in L. At energies where 
the angular momentum associated with the bounce, Tsr (A), exceeds LyiE) far fewer stars are captured by the SBH 
than in simulations that neglect the effects of GR. The iMerritt et al.l (I2011D study revealed that orbits experiencing 
the “bounce” were of such high eccentricity that their GR precession times were short compared with those of typical 
(i .e., less eccentric) stars at the same a. 

iHamers. Portegies Zwart fc MerritH (l2014f l coined the term “anomalous relaxation” to describe the behavior of orbits 
in this high-eccentricity regime, L < Lsb{E). Those authors verified the existence of the SB via an independent set of 
A^-body integrations, and also carried out test-particle integrations, using a much larger number of stars, from which 
they numerically evaluated the rates of diffusion in the anomalous regime. 

Based on these, and other, studies, two analytic expressions have been proposed for the location of the SB. The first 
compares the GR precession time with the time for the '/N torques to change L (jMerritt et al.l[2(TTTI ): 


’^sb(®) 



' M, Y 


N{a). 


(17) 


The second (iHamers. Portegies Zwart fc Merritt 12014 iBar-Or fc Alexander! 1201^ compares the GR precession time 
with the coherence time: 


^SB(a)^4- 


^coh(l^) 


(18) 


In spite of their disparate functional forms, the two expressions can yield numerically similar relations for 77.sB(a)> as 
illustrated below. The former relation appears to more accurately reproduce the barrier location in numerical studies 
to date; while the latter relation arises naturally w hen matching diffusion coefficients in the resonant and anomalous 
regimes (jHamers. Portegies Zwart fc Merrittil2014ll . 

Evaluating the former expression in the case of an unmodified Bahcall-Wolf cusp, n{r) oc yields 



where r™ is the radius containing a mass in stars of 2M,. The barrier as given by Equation (IT^ extends between the 
radii amin and Omaxi where 


^min 



^max 




(20a) 

(20b) 


The first relation follows from setting e = 0. The second relation is the intersection of Equation (fTOl) with the curve 
a(l — e) = Org, the periapsis of an orbit that intersects the loss sphere of radius ric = Org. Taking parameter values 
appropriate for the Milky Way: 

M, = 4 X 10® M©, 


TO* = 1, 0 = 15 
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yields 


amin ~ 6 — mpc, 

VPC/ 


^max 


140 



( 21 ) 


“Anomalous relaxation” is defined as angular-momentum diffusion of orbits with L < Lq'q^E). Paper I presented a 
derivation, based on a simple Hamiltonian model, of the diffusion coefficients in the anomalous regime: 


where T(a) = tcoh(a)/(A^)^ and 


(A7^)AR = ((A7^)2)AR = 

T T 

A = 1 M^,{a) a 

~ M, Tg' 


( 22 ) 


(23) 


The rapid, power-law drop predicted in the diffusion rates for TZ < TZsb is due to the adiabatic invariance of L under 
the effects of rapid precession. 

The derivation leading to equations (l^ - (l^ w as very approximate and one would l ike to verify those functional 
forms by comparison with A-body integrations. [Hamers. Portegies Zwart fc Merritt! (|2014f) attempted to do this. 
However, it was found that the value of N accessible to high-accuracy simulations {N < 100) was so small that the 
effects of anomalous relaxation could not be cleanly differentiated from the effects of classical relaxation at small L. 
Application of a more approximate, test-particle approach allowed Hamers et al. to increase the effective value of N 
by two orders of magnitude. The diffusion coefficients extracted from these experiments were found to be reasonably 
well described by equations (HHl-dlSl). 

In the present treatment, we account for the effects of anomalous relaxation by modifying the angular momentum 
diffusion coefficients m- We consider two sorts of modification with different functional forms: a power-law mod¬ 
ification, which reproduces equations (I22[) at small 72.; and an exponential modification, which implies a much more 
rapid decrease in the diffusion rate toward small 72. 


3.1. Power-law modification 

To account for anomalous relaxation, the angular-momentum diffusion coefficients of Equation (1131) are modified as 
follows: 


(A72) = (A72)ck +5i(f,7^)(A72)RR, {{AUf) = ((A72)2)ck + <?2(^,7^)((A72)2)rr. (24) 

The functions 51 (72) and ^2(72) should have certain properties. Both gi and 52 should tend to one as 72 —>■ 1. Since 

(A72)rr ^ 2A{£), ((A72)2)rr ^ 4A(f )72 

as 72 —>■ 0 f Equation I14|l . we require gi —>■ 72^ and <72 —t 72^ for small 72 so that the small-72 behavior of Equation 
(1^ is reproduced. The transition between the two regimes should occur at 72 ^ 72sb for both functions. 

An ad hoc functional form for <72 that satisfies these requirements is 


g2{£,'n) 



-2/n 


(25) 


where 722 ~ 72 sb('?)- The parameter n determines the rapidity of transition between the large-72 and small-72 regimes. 

The same functional form might be adopted for gi. Rather than make that choice, we first consider another possible 
constraint on 771 and ( 72 . 

The 72-directed flux coefficients that appear in the Fokker-Planck equation are given by Equations m- 


Dn = -{An) - ^{ASATZ) + ^-^{ASATZ) + ^^((A72)2) (26a) 

?«-(A72) + i^((A72)2), (26b) 

Dnn = l{{An)^). (26c) 

Since the 72-directed flux is 

df d f df 

~ ~ ^ ~^^^dTZ ~ 

it is reasonable to require that the diffusion coefficients in 72 satisfy 

Dn —t 0 , Dtztz —>• 0 (28) 
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at the boundaries TZ = {0,1}; in other words, that 

{m = ((A7^n = o 


(29a) 


at 7^ = {0,1}. Both the classical (Cohn-Kulsrud), and the resonant diffusion coefficients adopted here and in Papers 
I and II satisfy these conditions. 

Suppose that a stronger condition is imposed: D^i = 0, 0 < 7^ < 1. In this case, the Fokker-Planck equation 
describing diffusion in angular momentum reduces to 


dt dn dTZ\ ^^dn) 


lA 

2m 


(((ar)^)A) 


(30) 


which has a steady-state (zero-flux) solution f{TZ) = constant regardless of the functional form of {{ATZ)^). It could 
be argued that filZ) = constant is a reasonable form for a time-independent /, since it corresponds to an isotropic, 
or “maximum entropy,” state. The resonant diffusion coefficients adopted in Papers I and II satisfy this stronger 
condition. Setting D-jz = 0 also implies zero drift, i.e., zero flux in TZ in the absence of a gradient. For this reason, 
D-ji = 0 will henceforth be called a “zero-drift” condition. 

Returning now to the anomalous diffusion coefficients, we ask: what functional form for gi is required for zero drift? 
That is: 

Q = Dn = -5l(7^)(A7^)RR + i A [g^{TZ){{ATZf)nn] ■ (31) 

Assuming that 52 is given by Equation (1^51) . the result is 


9 i{TZ) = \i + 


R2{£) 


n 


-2jn 




2(l-7^) /7^2 


1 - 27^ 


n 


1 -b 


7?2(f) 


-in'! —2/n —1 


n 


(32) 


Since the zero-drift argument is one of plausibility only, we will consider a slightly more general expression for gi that 
includes “zero drift” as a special case: 


9im = 



2(l-7^)/7^l^ 


1 -b 


' Rijsy 

TZ 


-2/n-l 


(33) 


The only difference between Equations (1321) and (1331) is the introduction of a second parameter, TZi, in place of 7^2- 
This generalization still satisfies the zero-flux condition (^5)) at TZ = 0, regardless oiTZijTZ^- ki TZ = 1, gi and 92 
are both very close to one (especially since n will be chosen to be large) so that the resonant diffusion coefficients are 
recovered and the zero-flux condition is satisfied, again for any choice of TZi/TZ 2 . 

Correspondence of these expressions with the diffusion coefficients of Equation (1221) at small TZ would require 

7^?(£) = (6/5)rA(f), TZli£)=TAi£), (34) 


where 

rA{S)=i»iy-fy (Jfy. (35) 

To within factors of order unity, these relations imply 7^i « 7?.2 ~ R-sb (equation [TH]) . 

It is shown in the Appendix that at small TZ, these choices for gi and g 2 imply 

/ 7 ? \ ^ 7?2 

Dn^ 6 XA{£)i^—j , A=I-^ (36) 

so that the direction of the drift is determined by the relative sizes of TZi and TZ 2 , as follows: 

7^1 < 7^2 —^ TD'yi O 0 —y (j)ji > 0 (37a) 

TZi > TZ 2 —^ TDji > 0 —>■ (jjfi < 0 (3fb) 

where the expressions for 0^ assume f{TZ) = constant. Furthermore both the form of the steady-state f{TZ), and the 
steady-state flux (assuming the presence of a sink, i.e. that f{TZ) = 0 at 77. < TZq), depend sensitively on A, as shown 
in Figures ITOl and fTTl from the Appendix. In the zero-drift (A = 0) case, the steady-state flux is reduced by a factor 


r] Ki 2 



(38) 


compared with the flux that would obtain in the absence of anomalous relaxation (note that an empty loss cone has 
been assumed). For nonzero A, the reduction factor can be greater or smaller than this, tending toward a maximum 
value of one for large TZ 1 /TZ 2 (Figure fTTl). 
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3.2. Exponential modification 

IBar-Or fc Alexanded (|2014r) suggested a different functional form for ((ATT.)^) in the anomalous regime: 


((ATT.) ) oc exp- 


1 n. 


SB 


TT Ti? 


(39) 


an exponential cut-off toward small TT. While there does not seem to be strong support for this functional form in 
any published numerical simulations, we consider it here for completeness, and because it serves to illustrate how 
sensitively the evolution of / depends on the form assumed for the anomalous diffusion coefficients. 

Proceeding as above, we write 


Suppose 


(ATT) = (ATT)ck + hi(f,TT)(ATT)RR, ((ATT)^) = {{Mlf )ck + h2 (£, TT)((ATT)2)rr. 


h 2 {S,n) = exp ( 


(40) 

(41) 


where ~ TTsb(^’)- Since TT 4 <C 1, h 2 {S,l) ~ 1. The zero-drift condition would imply 


(ATT) = 2Ai£) 


1-2TT + 2(1-TT) ( ^ 


exp - 




We again generalize this expression by defining a second parameter, TT 3 , and writing 


hi{n) = 


1-b 


2(1-TT) /TT 3 


1 - 2TT 

At small TT, these expressions imply a flux coefficient 


TT 


exp - 


3. 


D 


n 


AAi£) 


TUV 

n) 


exp - 




l-|fexp|- 


Rl - TT2 ^ 1 
TT2 


(42) 


(43) 


showing that in this case, as in the power-law case, the sign of D-jz is determined by the relative sizes of TT 3 and TT 4 . 
Figures ITUl and [TT] illustrate the properties of the steady-state solutions /(TT). The behavior of / at small TT depends 
very sensitively on TT 3 /TT 4 . The reduction in the steady-state flux, for TT 3 = TT 4 , is 


T ]; 


o’^4 , 


TT 4 

TTo 




(44) 


3.3. Constraining the functional forms of the diffusion coefficients in the anomalous regime 

Based on the analysis presented above and in the Appendix, many properties of the steady-state solutions are 
expected to depend sensitively on the functional forms of the diffusion coefhcients in the anomalous regime, and in 
particular on the ratios TT 1 /TT 2 or TT 3 /TT 4 . 

In Paper I, the functional forms of the di ffusion coefficients in the resonant r e gime were successfully constrained 
by comparison with the numerical results of [Hamers. Portegies Zwart fc MerrittI ( 20141) . Those authors used a test- 
particle integrator to extract values of the angular momentum diffusion coefhcients for stars orbiting near a SBH in 
nuclei with n(r) oc and r~^ . 

The Hamers et al. integrations included post-Newtonian terms in the equations of motion, both for the held and test 
stars, and the suppression of angular momentum diffusion below the Schwarzschild barrier was clearly seen. Figure [I] 
provides an illustration: it shows the hrst- and second-order diffusion coefhcients for stars in a single energy bin, 

in integrations of a model with n(r) oc r~^ EI_ Overplotted are analytic diff usion coe fhcients from the two families 

considered above: power-law (Equations lAll IA2|) and exponential (Equations lAll IA23I) . 

These hgures, and similar ones made for stars at other energies, motivate the following conclusions (some of which 
were presented already by Hamers et ah): 

1. An e xponential dependence of the diffusion coefhcients on TT in the anomalous regime (|Bar-Or fc Alexandeii 
1201411 is ruled out. 

2. The power-law dependence of (ATT) and ((ATT)^) on TT predicted in Paper I, and reproduced here in Equations 
(EH), is conhrmed, particularly in the case of the second-order coefficient for which the noise is smallest. 

3. The value of TT that dehnes the transition between the resonant and anomalous regimes, called here TT 2 , is well 
predicted by Equation (fTSl) . 

4 These data were kindly provided by A. Hamers. 



























R R 

Fig. 1.— Fits to diffusion coefficients extracted from the test-particle integrations of lHamcrs. Portcgics Zwart fc McrrittI 1)20141') . for the 
radial bin (a) = 11.9 mpc. In the upper panels, red symbols are — {A7^}. Diamond symbols a re b i nned data, with errors, to 
guide the eye; fits were based on the unbinned data. Left: analytic curves are Equations iEH i. iix^ . the power-law model for 
anomalous relaxation. 7^2 = 0.110 (shown by the dashed line i n the b ottom panel); in the upper panel, 7^i/7^2 = {1,1.5, 0.7}. Right: fits 
of the exponential model for anomalous relaxation. Equations llAll l. 1IA23II . 7^4 = 0.080, 7^3/7 ^4 = {1,1.25, 0.8}. Vertical dotted lines are 
estimates of the value of TZ at which classical relaxation dominates anomalous relaxation (Eos. fSTll . assuming that the power-law forms of 
the anomalous diffusion coefficients are correct; the analytic for ms would not be expected to describe the data below these lines. In the 
lower-right panel, the additional dotted line at 7^ 0.33 is Eq. Esji, an estimate of where classical relaxation would begin to dominate if 

the exponential forms of the anomalous diffusion coefficients were correct. 


An attempt was made to find the best-fit value(s) of 7^i/7^2 by searching for parameter values that optimized the 
fits to data like those in Figured] Unfortunately, the results so obtained were found not to be robust. This was due in 
part to the greater noise associated with data below the Schwarzschild barrier, particularly in the case of the first-order 
coefficient. In addition, the much greater variation in the amplitudes of (A7^) and {{ATZ)'^) at a single energy meant 
that the best-fit solution depended sensitively on the relative weighting of the data at different values of TZ. In the end, 
all that could be concluded was that the first-order diffusion coefficients are consistent with TZi —TZ 2 , the “zero-drift” 
condition; but that values of TZi /7^2 moderately greater or less than one could not be excluded using these data. 

The fact that the steady-state form of f{TZ), and the loss-cone flux, depend sensitively on 72.i/7?.2 suggests a sec¬ 
ond way to constrain the anomalous diffusion coefficients: insert them into the Fokker-Planck equation and integrate 
forward from initial conditions like the ones used by Merritt et al. (2011) in their small-iV simulations. As discussed 
in Hamers et al., the value of N in those simulations was too small to allow direct extraction of the diffusion coeffi¬ 
cients. However, the time-averaged, or integrated, properties of the V-body models were reasonably well determined, 
particularly given that multiple 8) realizations of the same initial conditions were integrated, allowing the variance 
in the results to be reduced by averaging. 

The IMerritt et all (120 1 lH initial conditions consisted of 50 stars, of mass m* = 50Mq, distributed as n(r) oc r~^ 
around a SBH of mass M, = IO^Mq. The initial distribution was truncated for orbits with semi-major axes above 10 
mpc and below 0.1 mpc. Integrations were carried out both with, and without, post-Newtonian terms in the equations 
of motion, up to order 2.5 PN. Capture of stars by the SBH was allowed to occur whenever the orbital periapsis fell 
below 8rg = 8GM,/c^. Each realization of the initial conditions was integrated for a time corresponding to 2 x 10® yr 
(with PN terms) and 10^ yr (without PN terms). The average capture rate in the relativistic integrations was about 
one event per 10® yr; in the absence of the PN terms, mean capture rates were about a factor 20 higher. 

There is no ambiguity in representing the Merritt et al. iV-body initial conditions as a smooth f{£,TZ). However, 
the Fokker-Planck algorithm has a number of parameters that must be specified, in addition to those that define the 
anomalous diffusion coefficients (n and TZi/TZ 2 , or 'IZz/TZ^j. Those parameters include 

^min; In A, 


At, Nx, Nz- 


(45) 
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fmin is the binding energy at the edge of the {£,TZ) grid; it should be small enough that few stars diffuse to £ < Smin 
during the course of the integration. The value 10“®c^ was chosen, which is the energy of an orbit with semimajor 
axis a = 5 X lO^rg « 2.5 pc, or ~ 250 times the maximum a-value of the initial conditions. The number of grid points 
was Nx = Nz = tS'i. The quantity In A was set to 15 in most of the integrations, except for one set in which smaller 
and larger values (from 11 to 19) were tried. The Coulomb logarithm only appears in the expressions for the classical 
diffusion coefficients; since evolution of these models is dominated by resonant relaxation, the results are expected to 
be weakly dependent on In A and this was found to be the case. The integration time step. At, was set to 2000 yr, i.e. 
10^ steps per integration. 

A natural choice for the parameter rxcjrg in the Fokker-Planck integrations would be 8, the same value assumed in the 
A-body integrations. In the case of “plunges” - captures that occur without significant energy loss due to gravitational 
radiation - this would be the correct choice. However, some of the A-body capture events were “EMRls,” for which 
capture was preceded by significant energy loss due to the 2.5PN terms. No such loss terms are included in the 
Fokker-Planck integrations described here. Roughly speaking, the effect of the 2.5PN terms is to shift the location 
of the loss cone toward larger 7^ (i.e. larger orbital periapsis) at each £ (see e.g. Figure 5 of Merritt et al. 2011). 
To evaluate the effect on the Fokker-Planck results of ignoring the 2.5PN terms, a set of integrations was carried out 
setting Hc/'Cg = 32, four times its value in the A-body integrations. 



time (yr) time (yr) 

Fig. 2.— Time-avera ged loss rates , defin ed as the total number of stars lost until time divided by t. (Blue) squares are fro m the 
A^-body integrations of fMerritt et all II201U) (Figure 2c of that paper); curves are from the Fokker-Planck integrations described in gM] 
The sudden jumps in the A^-body loss rates at early times correspond to single capture events. Left (right) panel shows 
results using the power-law (exponential) forms of the anomalous diffusion coefficients in the Fokker-Planck code. Left: the six sets of 
curves are for 7^i/7^2 = {2,1.2,1, 0.9, 0.8, 0.7}, from top to bottom. The different line styles are explained in the text. Right: the six 
sets of curves are for = {1.4,1.2,1.1,1, 0.95, 0.9, 0.7}, from top to bottom. Curves for '7Zil'7Z2 = 1 or 7^3/7^4 = 1 (“zero drift”) are 

thicker. 


Figure [2] compares time-averaged loss rates in the V-body and Fokker-Planck models, defined as the total number of 
stars lost until time t d ivide d by t . The left panel shows results using the power-law forms of the anomalous diffusion 
coeffi c ients, Equations iAil), (|A2I) . with n = 8; the right panel shows results using the exponential forms, Equations 
(|A1I) . (IA23|) . Each panel shows results for several values of TZi/TZ 2 (power-law) or TZ^/TZi (exponential), as specified in 
the caption. For each value of this ratio, three integrations were carried out, varying the way in which 72.2 or 72.4 were 
related to 72sb- One of the three integrations (shown by the solid curves) equated 722 or 724 with TZg^ , Equation (fl^ . 
The other two integrations adopted larger or smaller values: by a factor two or one-half (in the power-law models), or 
by a factor or 1/v^ (exponential). These additional integrations are shown by the dash-dotted curves in Figured 
Larger (smaller) values of 722 or 724 generally resulted in smaller (larger) loss rates, at least at early times. 

Figure [2] suggests that both functional forms of the anomalous diffusion coefficients are able to reproduce the iV-body 
capture rates, as long as 72i/722 or TZ^lTZ^ are not too different from one. The best correspondence is achieved, in 
both cases, when this ratio is slightly less than one, and this result remains unchanged even when the values of 722(f) 
or 724 (f) are substantially modihed. Recall that 72i < 722 or 723 < ’^4 imply D-ji < 0, i.e. (jin > 0, i.e. a drift toward 
larger 72 1 Equation 15711. 

Changing the parameter n in the power-law diffusion coefficients from n = 8 to n = 32 had almost no discernible 
effect on the loss rates. Varying In A or ric/rg in the amounts described above did result in noticeable changes, but by 
amounts comparable with the ranges shown in Figure!^ due to variations in the definition of 722 or 724. In every set 
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of integrations, correspondence with the 7V-body loss rates was best for TZi/TZ 2 < 1 


or 72 . 3 / 72.4 < 1. 



■ 

; 

/ / / / / JP \ 

11 / /° / p 4 
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ogio(^-e) 
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/ / / /, 
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= logio(1-e) 


X = logio(1-e) 


X = logio(1-e) 


Fig. 3.— Time-avera ged angular momentum distributions for stars at a single energy. The symbols are from the A^-body integrations 
of IMerritt et al.l d201U') (Figure 11 from that paper); the triangles exclude stars that eventually became EMRIs while the squares include 
those stars. Curves were extracted from Fokker-Planck integrations that used the power-law (top) or exponential (bottom) forms of the 
anomalous diffusion coefficients; the curves show where £^4 is the energy corresponding to orbits of semimajor axis 4 mpc. The 

A^-body data were computed using stars with instantaneous a values in a range Alog^pU = ±0.05 centered on a = 4 mpc. Curves are 
distinguished by the value of 7 ^i/ 7^2 = {0.7, 0.8, 0.9,1,1.2,1.5, 2} (top) or 7 ^ 3 / 7^4 = {0.7, 0.9, 0.95,1,1.1,1.2,1.4} (bottom); the larger this 

ratio, the larger the value of / at small X. The middle panel used 7^2 = ~ right panels used 7^2 = 

and 7^2 = ^ ^ respectively (top) or 7^4 = "^4 — ^SB^ ^ (bottom) Filled circles mark 7^ = 7^2 (top) or 7^ = 7^4 

(bottom). 


Figure [3] makes another comparison between Fokker-Planck and A^-body models. Plotted there are time-averaged 
angular momentum distributions at a single energy. These are displayed as dN/dX, where X is defined as in Figure 
11 Of IMerritt, et, all (IMIh : 

X = logio (1 - e) (46) 

with e = \/l — 72 the orbital eccentricity. Figure [3] implies that values of 72i/722 or 723/724 of unity (“zero drift”) 
or greater can be securely ruled out - consistent with Figure [2l In the case of the power-law forms of the diffusion 
coefficients (upper panels), the best correspondence with the V-body results seems to occur for 

722 ~ 272^“\ 0.8<72i/722 < 0.9. (47) 

In the case of the exponential forms of the diffusion coefficients (lower panels), correspondence with the V-body results 
seems to require 

722 ~V^72^b^ 0.9 <72i/722 < 0.95. (48) 

Once again, the best correspondence is achieved with diffusion coefficients that imply a non-zero drift, in the direction 
of increasing 72, in the anomalous regime. 

Figure m shows representative plots of f{TZ;£) at one £ from the Fokker-Planck models at the final time step (roughly 
2.5 X 10® yr). Dashed curves show models having parameters similar to those found to correspond best to the 7V-body 
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Fig. 4.— Plots of /(T?.) at the final time step of the Fokker-Planck integrations described in i)3.3l for stars at a single energy, corresponding 
to orbits with semimajor axes of 4 mpc. The left panel shows integrations based on the power-law form of the anomalous diffusion 

(a) 

coefficients; the right panel is based on the exponential form. On the left, the parameter 1^2 = , and the different curves correspond 

to 'R.iI'R.2 = {0.7, 0.8, 0.9,1,1.2,1.5, 2}; the larger this ratio, the larger the value of / at small TZ. On the right, Tti = \/n'R.g^ and 
Tiil'R-A = {0.7, 0.9, 0.95,1,1.1,1.2,1.4}. The thick solid curves in both panels are for 'R,\JTi .2 = 1. The dashed curves show models that 
were judged to best reproduce the W^-body data, based on results like those in Figures [2] and |3] Filled circles mark 'R. = 7?.sb; verti cal 
dott ed lines show the values of Ti. below which classical relaxation is expected to dominate anomalous relaxation at this energy (Eqs. 1511 
and l53l l. 


results. These solutions always exhibit a rapid drop in the steady-state f{TZ) below the Schwarzschild barrier. That 
drop would be even steeper in the absence of classical relaxation, which dominates the evolution in TZ at small TZ (the 
region below the vertical dotted lines in the figure), thus maintaining a relatively high diffusion rate at small TZ. 

One final argument can be made in support of diffusion coefficients that satisfy D-j^ <0. As described in 
lAntonini fc Merritt! (|2013ll . stars orbiting near the Schwarzschild barrier are observed to exhibit a “buoyancy” phe¬ 
nomenon: should they cross the barrier from below (TZ < TZsb) to above {TZ > TZsb), they tend to remain above, and 
vice-versa. This behavior is consistent with a drift toward larger TZ, as implied by D-jz < 0. 


3.4. Dominance of classical relaxation at small TZ 

In the regime of anomalous relaxation {TZ < TZsb{£))■, timescales for angular momentum diffusion become long for 
small TZ. One consequence is that classical relaxatio n, which (by assumption) is not affected b y precession, can once 
again dominate the evolution in angular momentum (|Hamers. Porteeies Zwart fc Merrittll20f^ . 

We estimate the value of TZ at which this occurs by equating the first and second terms on the right hand sides of 
Equation (1241) . We simplify the expressions by using the limiting forms o f the diffusion coefficients as 7^ —» 0. _ 

For the classical diffusion coefficients. Equations (11) and (12) from [Hamers. Portegies Zwart fc MerrittI (j2014fl . 
together with the transformation Equations (32) from Paper I, yield 


{ATZ)ck^ {£), 


A^{a) 


In A 

Cnrr(7) 


{{ATZf) CK —>■ 2TZA^{£), 
N{a) 

P{a) 


(49a) 

(49b) 


where the symbol denotes the limit of small TZ and the function C'nrr(7) is given in Appendix B of Hamers et 
al. (2014); their calculation assumes p{r) oc r~^. 

For the anomalous diffusion coefficients, Equations (l24l) . (|2^ and (1151) give the limiting forms 


{ATZ)ab.^6A{£) 



((A7^)2)AR ^ 4A{£)TZ 



(50) 


in the power-law case, with A{a) defined in Equation (fT^ . Equating (A7?.)ck with (A72.)ar or ((A72.)^)ck with 
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((A7e)2)AR yields 


_ InA P 

6 CnrO^ teoh 

7^2 _ InA P 
T^l 2 Cnro;^ tcoh 


(1st order), 


(2nd order). 


Replacing Pi and P 2 by 


in Eqs. and setting a = 1.6 yields 


7^2 = (1.0,3.1) 


InA tcoh 

Cnr V a / P 


(51a) 

(51b) 


(52) 


where the constants in pare ntheses refer to first- and second-ord e r diffu sion coefficients respectively. Equation (I52|) is 
similar to Equation (23a) of lHamers. Portegies Zwart fc Merritt) (120141 ). 

Equations ((CTl were used to plot the vertical dotted lines in the left-hand panels of Figures [Hand |4] (note that these 
two figures refer to different mass models). In Figure |TJ the line predicts reasonably well the value of P at which the 
data begin to deviate from the fitted curves. In the case of Figure SI the effects of classical relaxation can be seen in 
the steady-state f{P), which drops more gradually to zero below the SB than it would if only anomalous relaxation 
were acting (Figure [TO]). 

If the exponential forms of the anomalous diffusion coefficients are correct, then the value of P at which ((A7 ^)^)ck = 
((A7^)2)AR becomes 


P^ 


log 


2oi (7 nR ^coh 

InA 


(53) 


This value for P is plotted, in addition to the value given by Equation (l5^ . on the lower right-hand panel of Figure |TJ 
Note that under this hypothesis, the range in P over which anomalous relaxation would be relevant would become 
very small; in effect, classical relaxation would dominate the evolution everywhere below (and sometimes even above) 
the Schwarzschild barrier. We reiterate that there is no support for the exponential form of the anomalous diffusion 
coefficients in any published numerical simulations. 


4. STEADY-STATE SOLUTIONS 

Paper II presented steady-state solutions for f{£,P) obtained via integrations of the Fokker-Planck equation with 
diffusion coefficients given by Equation (USD (no anomalous relaxation) and outer boundary condition 

r{£P^,P\e) = r{£P^,P\Q) (54) 

with / the phase-space density and £min the minimum value of £ on the energy grid. (Asterisks denote dimensionless 

quantities; see Eauation ll2l and Paper I.) Initial conditions for / were based on an isotropic power-law model, n oc , 

f oc f but with a simple modihcation to account for the presence of the loss cone: 

/(£:,7^,^ = 0 ) = 0 , P<Pic(£). (55) 

Among the parameters that were varied in the integrations of Paper II were m*/M, and the initial density at large 
radii; the latter was chosen to have one of three values, bracketing the estimated value for the nucleus of the Milky 
Way. The physical radius of the loss sphere, ric, was hxed at 15rg = 15G'M,/c^, roughly the tidal disruption radius of 
a solar-type star. 

We repeated a subset of those integrations, now using the modified expressions for the angular-momentum diffusion 
coefficients that account for anomalous relaxation: either power-law ('Equation 1^^ or exponential (Equation 0DD . 

The dimensionless parameter was set to 2.5 x 10“^. Assuming M, = 4 x IO^Mq, the (dimensional) stellar 

mass becomes m* = I.OM 0 . The outer boundary condition was chosen, as in Paper II, to give one of the following 
three values for the mass density at one parsec: 

{1.9 X 10"^, 3.5 X 10®, 6.1 X 10®}Mqpc“® (56) 

where again M, = 4 x IO^Mq has been assumed. For each of these parameter choices, two other parameters that 
appear in the expressions for the anomalous diffusion coefficients were varied: 

1 . P 1 /P 2 (power-law) or P^/Pi (exponential) ; 

2. The ratio P 2 /PSB or Pi/Psb- 

Based on the results described in the previous section, the following parameter values were considered: 

• 7 ^l/ 7^2 = ( 0 . 8 , 1 . 0 , 1 . 2 }; P 3 /P 4 = (0.9, 1 . 0 , 1 . 1 } 

• 7^2/7^sB = (0.5,1.0,2.0}; Pa/Psb = (0.5,1.0,2.0} . 
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In R 



In R 



In R 


Fig. 5. — Angular momentum diffusion coefficients at t = 0 in the models of [j4] Top: {A7?.}ck (left), {(A7?.)^}ck (right)- Middle: 
{A7?.)ck + (^I^)rr (left), ((A7^ )^)c k + ((^1^)^)rr (right). Bottom: Modified forms of the diffusion coefficients that account for 
anomalous relaxation. Equations II24II . with parameters as given in the text. The physical loss cone radius is shown as the thick (blue) 
curve in each panel; / = 0 is assumed inside this curve at t = 0 (“empty loss cone”). The thin (red) curve that lies inside the loss cone 
is the quantity TZf)(S) defined in Paper I. In the panels showing (AT?.), the red contours indicate —(AT?), i.e. (AT?) < 0 in these regions. 

In the lower panels, two expressions for the location of the Schwarzschild barrier, 'Rg^{£) and T?.!*^ (£), are shown respectively as the thin 
and thick magenta curves. The dashed curves in the lower panels indicate where the timescales for classical and anomalous relaxation are 
equal. Contour values are the same in all frames. 


For 72 .sb(^’), the expression (11811 was used. The parameter n that appears in the power-law form of the anomalous 
diffusion coefficients was set to 8 in all integrations. 

Figure [S] plots the diffusion coefficients {ATZ), {{ATZ)"^) as computed by the code at t = 0, in models having the 
middle of the three values for the mass density at one parsec (Equation [35]). The top two frames plot the classical 
diffusion coefficients: 

(A7e)cK, ((A7^)')CK. 
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The middle two frames plot diffusion coefficients that account for resonant relaxation: 

(A7^)ck + (A7?.)rr, ((A7?.)^)ck + ((A7^)^)rr. 

These are the same expressions adopted in the integrations of Paper 1. The lower set of frames show the diffusion 
coefficients of Equation (1241) , which account for anomalous relaxation (power-law modification, with n = 8, 0-2 = 27^gg , 
and 7^1 = 7 ^ 2 )- In the case of the exponential modification (not shown in this figure), the diffusion coefficients drop 
off more steeply below the Schwarzschild barrier. However, this drop is mediated, in all models, by the presence of 
classical diffusion, which dominates again at sufficiently small 72., as discussed above. 



n In i? In i? 

Fig. 6. — Grey-scale of log/* in a set of six, steady-state models. Top: models integrated using the power-law forms of the anomalous 
diffusion coefficients, with 7^2/’7 ^sb = 2, and with 7^i/7^2 = 0.8 (left), 1.0 (middle) and 1.2 (right). Bottom: models integrated using the 
exponential forms of the anomalous diffusion coefficients, with 7^4 /'^sb = 2, and with 7^3/7^4 = 0.9 (left), 1.0 (middle) and 1.1 (right). 
Other parameters are given in the text. Curves have the same meaning as in Figure (5] 


Steady-state solutions, /*(£■*,72*), are shown in Figure [H again for models with p(lpc) ~ 3.5 x lO®M 0 pc“^. The 
top(bottom) panels show solutions obtained using the power-law(exponential) expressions for the anomalous diffusion 
coefficients, with 722/72 sb = 2 or 724 / 72 sb = 2. What varies, from left to right, is the choice of 72i/722 (top) or 723/724 
(bottom). When these ratios are unity (“zero-drift”), the steady-state solutions are characterized by /(72) ^ constant 
near the SB. When these ratios are greater or less than one, the steady-state solutions behave in roughly the way seen 
in Figures l4l and flOl becoming either greater or smaller in the region below the Schwarzschild barrier, before dropping 
to zero at the loss cone boundary. Evidently, the form of the steady-state / in this region depends very sensitively on 
deviations of that ratio from unity. Depending on the value of that ratio, / in the region below the barrier can either 
be strongly depleted, or strongly enhanced, compared with the “zero-drift” solution. 

Angular-momentum-averaged distribution functions, defined as 

7{£)= [ f{£,n)dn, ( 57 ) 

JTliaiS) 

are plotted in Figure [7] for all of the steady-state solutions. Shown for comparison, as the thick solid curves, are / 
for models computed without anomalous diffusion; these are the same three curves plotted in Figure 5 of Paper II. 
As discussed in that paper, inclusion of the resonant diffusion coefficients has the effect of sharply truncating the 
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Fig. 7. — Angular-momentum-averaged distribution functions, f{S), for all of the steady-state models; left(right) panel used the po wer - 
law(exponential) forms of the anomalous diffusion coefficients. The three values of the large-radius density normalization, Equation IMJ, 
are indicated by the three colors (high density, magenta; intermediate density, black; low density, blue). In each of the three sets of models, 
the solid curves have 'R,i/'R ,2 = 1 (left) or TZilTZi = 1 (right); the dashed curves have 'R.xj'R.^ = 1-2 (left) or = l.f (right); the 

dotted curves have 'R,i/'R ,2 = 0.8 (left) or = 0.9 (right). For each choice of these parameters, there are three curves, with the same 

line style, corresponding to the three choices {0.5,1.0.2.0} for 7?.2/7 ?.sb or 7?.4/7 ?.sb. The three thick, solid curves in each panel are from 
integrations that did not account for anomalous relaxation; these are the same curves plotted in Figure 5 of Paper II. 


steady-state f{£), at binding energies above a certain value, where the timescale for resonant diffusion (in angular 
momentum) drops below the timescale for classical diffusion (in energy), and stars are carried rapidly into the SBH. 
The truncation of / can still be seen in the new models; but it is mediated by the presence of anomalous relaxation. 
The reason is the increase in the angular-momentum diffusion time when anomalous relaxation is accounted for: stars 
“pile up” near the Schwarzschild barrier, until their density is high enough to drive the requisite flux. 

Closely related to f{£) is p(r), the mass density. Figure |8] shows steady-state density profiles for all the integra¬ 
tions. Also shown are three density profiles from Figure 4 of Paper II (no anomalous relaxation), which exhibit cores 
corresponding to the depletion in f(£) at large binding energies due to resonant relaxation. Once again, the lesser 
depletion in the models that account for anomalous relaxation translates into cores of lesser prominence. Indeed in 
the models with TZi < 1^2 or TZz < 72.4, the steady-state density profiles turn out to be very close to the classical 
Bahcall-Wolf cusp at all radii plotted. This plot confirms a conjecture made in Paper I: namely: that the inhibition in 
angular-momentum diffusion associated with the Schwarzschild barrier would reduce the ability of resonant relaxation 
to form a core. The generality of this result is discussed in ^ 

Even models having similar p{r) or f{£) can have very different loss rates into the SBH, since the latter depends 
also on the timescale for angular-momentum diffusion. The flux of stars into the loss cone, Equation (I^, is plotted 
for the steady-state models as a function of energy in Figure IHl Shown for comparison are loss rates in steady-state 
models computed using only classical, or classical plus resonant, relaxation. These plots show that the flux of stars 
into the SBH can depend strongly on the assumed forms of the anomalous diffusion coefficients. There are two, 
competing effects. Including anomalous relaxation tends to increase the steady-state density at small radii compared 
to p{r) computed using resonant relaxation alone, resulting in a larger flux. Anomalous relaxation also increases the 
timescales for angular momentum diffusion, which tends to reduce the flux. 

The dependence of the flux on the assumed value of 72i/722 or 72 . 3 / 72.4 is similarly complex. At low binding energies. 
Figure IHl shows that small values of this ratio imply lower fluxes; while at high binding energies, the reverse is true. 
The former result is consistent with the analysis in the Appendix (Figure fTTI) . The reason for the latter result can be 
seen in Figure 0 models with smaller values of this ratio maintain larger / at large binding energies, which tends to 
increase the flux. Total, or integrated, loss rates for these models can be computed using Equation ([8]) . 
The results turn out to be nearly the same within a few percent — for each of the models, roughly 
7 X 10“^ stars per year. As discussed in Paper I, this is a consequence of the fact that the total loss rate 
is dominated by stars at low binding energies, where the effects of resonant and anomalous relaxation 
are small. 
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radius (pc) 


radius (pc) 


Fig. 8. — Mass density as a function of radius for the steady-state models. 
Figure fn Scaling to physical units assumes M» = 4 x 10 ®Mq. 


Line colors and line styles have the same meanings as in 




Fig. 9. — Dimensionless flux of stars into the lo ss c one as a function of energy in the steady-state models; all of these models used the 
intermediate of the three values given in Equation ll56t for the density at one parsec. Left(right) panels adopted the power-law(exponential) 
forms of the anomalous diffusion coefficients. The two, thick solid curves in each panel show models that included only classical relaxation 
(curves that peak near the right) and only classical plus resonant relaxation (curves that peak near the left). Line styles have the same 
meanings as in Figures [7] and js] Circles are plotted at values of E corresponding to the outer edge of the Schwarzschild barrier. 


5. DISCUSSION 

5.1. Steady-state density profiles 

In Paper II, the formation of cores due to resonant relaxation was discussed. Equation (38) of that paper gave an 
estimate of the radius of the core so formed: 


0.028 


In A 

IF 


-4/5 


(58) 
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Here, is the gravitational influence radius of the SBH, defined as the radius of a sphere containing a mass in stars 
of 2 M,. 

As described in this paper, at least in the case of the Milky Way, the inclusion of anomalous relaxation tends to 
counteract core formation by causing stars to accumulate near the Schwarzschild barrier. 

Here we consider the generality of that result. Suppose that nuclei of galaxies fainter than the Milky Way have 
SBHs that satisfy the M, — a relation: 

M, _ ( a 

■ VlOOkms"^ 

(|Merrittll2013L Eq. 2.33) and that their nuclear densities are close to the Bahcall-Wolf form, p{r) oc . The latter 
is not too different from p oc r“^, for which « r/j = GM,/a^. Under these assumptions, 

pc, M. < 4 x 1 O 6 M 0 . (59) 

We first verify the assumptions that led to Equation (1551) . That equation was derived assuming that tcoh = tcoh.M 
(Equation [HI). The radius at which tcoh,M = tcoh.Sj for a nucleus with p{r) oc is given by Equation (30) from 

Paper H: 



^coh 




1.2 X 10"^ 


M, 

IO^Mq 



1.3 X 10"^ 


M, 

IO6M0 


0.18 


(60) 


suggesting that indeed fcoh Aoh.M at a = Ocore in the galaxies of interest. 

Next we ask how Ocore compares with the radii that define the Schwarzschild barrier. Equation (1201) gave approximate 
limits on the extent of the SB in a p oc nucleus, which we recast here as 
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Using Tg ~ 4.78 x 10 ® pc, these become 
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Setting Oinax < Ocore implies that anomalous relaxation is unlikely to affect the formation of the core due to resonant 
relaxation. This condition is: 
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(64) 


In the case of the Milky Way (M, si 4 x lO^M©), satisfying this condition for m* = M© would require > 40 pc - 
about ten times larger than the value inferred from stellar kinematics. This is consistent with Figure lU which showed 
that anomalous relaxation inhibits the formation of a core for all reasonable values of In the case of galaxies 
with central black holes less massive than the Milky Way’s, Equations (| 6 ^ and (l59ll allow us to write the condition 

(^max ^ ^core 

0.48 / \ 0.48 /I » \ -0.86 


M, < 5 X 10® 


_0y-4« /^y-48 


15 




V 15 / 


M«. 


(65) 


Thus, the core formed by resonant relaxation is expected to become progressively more prominent as M, is reduced be¬ 
low its value in the Milky Way. This fact is likely to be important in determining the rate of formation of gravitational- 
wave sources, particularly since theoretical estimates often focus on galaxies with M, < lO^M©. Estimating this rate 
will be the topic of upcoming papers in this series. 
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5.2. Constraining the forms of the anomalous diffusion coefficients 

Until recently, discussions of gravitational encounters near a SBH have usually been presented in terms of dif¬ 
fusion timescales', that is; in terms of second-order diffusion coefficients like ((AL)^) (iRauch fc Tremainel llQ^ _ 

iHopman fc AlexandeiJ2006allbl : iGiirkan fc HopmaiJ2007t lEilon et ^120091 : iMadigan et al.ll2011ll . iHamers. Portegies Zwart k, Mm 

(|2014f) were apparently the first to consider the forms of the first-order diffusion coefficients. Of course, both first- and 
second-order diffusion coefficients are essential when computing the evolution of / via the Fokker-Planck equation. 

As shown here through a number of examples, the form of the steady-state f{E,'R-) near the Schwarzschild barrier 
can depend very sensitively on the relative amplitude of the first- and second-order coefficients in the anomalous- 
relaxation regime. A natural case to consider is that of “zero drift”, in which the two diffusion coefficients imply a net 
flux in angular momentum that is zero when f{TV) = constant 1 ^3.II) . While it may be natural - it is consistent with a 
“maximum-entropy” steady state - this assumption is not compelled by any physical argument of which we are aware. 

In stellar dynamics, incorrect conclusions drawn from entropy arguments are legion, and numerical experiments, when 
available, would seem to be a bette r guide. As discuss ed in detail in 93.31 , the highest-quality A^-body simulations 
carried out to date of this regime (jMerritt et al.l[2f)TTh seem to require anomalous diffusion coefficients that differ 
slightly, though significantly, from the “zero-drift” condition. A state of “positive drift” seems to better characterize 
the existing simulations. 

Elucidation of the long-term effects of gravitational encounters in the Schwarzschild regime near a SBH will ultimately 
require a better specification of the anomalous diffusion coefficients. By far the best way to do this - at least in principle 
- is via direct A^-body integrations, which impose the fewest approximations. In practice, integrations of the required 
accuracy become very time-consuming when N > 10^. A major effort should be devoted to increasing the efficiency 
of the iV-body integrators. 


6. SUMMARY 

Integrations of the Fokker-Planck equation describing f{E,L,t), the phase-space density of stars around a super- 
massive black hole (SBH) at the center of a galaxy, were carried out using a numerical algorithm described in two 
earlier papers (|Merrittll2015alfO l. Diffusion coefficients describing classical, resonant and “anomalous” relaxation were 
included; the latter accounting for the evolution of orbits in the regime below the Schwarzschild barrier (SB) where the 
timescale for general relativistic precession is s hort compared with the coherence time, invalidating the assumptions 
that underlie the theory of resonant relaxation (iMerritt et al.l[2011I) . The principal results follow. 


1. Since a good theoretical understanding of anomalous relaxation is still lacking, two functional forms were considered 
for the angular momentum diffusion coefficients in this regime, having either a power-law or exponential dependence 
on TZ = Lf /Lf. In either case, a further choice must be made in terms of how to relate the first- and second-order 
diffusion coefficients. It was argued that a natural starting point is a “zero-drift” condition that implies no net flux 
in angular momentum when f{TZ) is constant. Parameterized functional forms for {A.TZ) and ((A7^)^) were proposed 
that have the “zero-drift” property as a special case. 

2. Two attempts were made to cons train the forms of the anomalous diffusion c oefficients by comparison with published 
numerical simulations. First, as in [Hamers. Portegies Zwart fc Merritt! (I2014f) . diffusion coefficients extracted from a 
large set of test-particle integrations were compared with the two functional forms. The power-law form was found 
to be strongly preferred, at least in the case of the second-order coefficient, confirming a result already presented in 
that paper. The first-order coefficient was also well fit by the power-law form, although data were more noisy and 
no clear preference could be established for the “zero-drift” conditions. Second, a set of Fokker-Planc k integrations 
were c arried out based on the same initial conditions that were used in the exact A-body integrations of iMerritt et al.l 
(|2011f) . Two properties of those models were then compared: the dependence of / on TZ, and the capture rate; in 
both cases, results from the A-body integrations were averaged over a set of different runs to reduce noise. Both the 
power-law and exponential forms for the diffusion coefficients could be made consistent with these data; it was argued 
that this was due in part to the effects of classical relaxation, which always dominates the diffusion rate at sufficiently 
small TZ. However, a clear preference was established for diffusion coefficients that imply a steady-state drift toward 
larger TZ, inconsistent with the “zero-drift” hypothesis. 

3. Fokker-Planck integrations were then carried out to find steady-state models having parameters similar to those of 
the nuclear star cluster in the Milky Way. These models were identical to the steady-state models computed in Paper 
H except for the inclusion of the anomalous diffusion coefficients. The steady-state f{£,TZ) in regions of phase space 
below the Schwarzschild barrier {TZ < TZsb) was found to be most strongly dependent on the assumed relation between 
first- and second-order diffusion coefficients. Diffusion coefficients satisfying the “zero-drift” condition produced steady- 
state solutions in which / was nearly constant with respect to TZ below the SB. Integrations incorporating positive- or 
negative-drift diffusion coefficients had steady-state /’s that respectively increased or decreased below the SB, before 
falling to zero at the loss cone boundary. In all of these models, departures of the steady-state density, n(r), from the 
classical Bahcall-Wolf solution were less pronounced than in the models of Paper H that did not incorporate anomalous 
relaxation (a result that was suggested already in that paper). The reason is the tendency of stars to accumulate near 
and below the SB, thus counteracting the depletion that occurs when only resonant relaxation is accounted for. 
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4. Steady-state loss rates in the presence of anomalous relaxation differ from those in all models previously published, 
for two reasons: different steady-state phase-space densities, and different forms of the angular-momentum diffusion 
coefficients. A robust conclusion is that the incorporation of anomalous relaxation implies a lower capture rate at 
energies where the SB exists, compared with models that only incorporate resonant relaxation, and this is true in 
spite of the generally higher steady-state densities in the former models. However the reduction in the capture rate 
depends sensitively on the parameters adopted for the anomalous diffusion coefficients, being most(least) extreme for 
the positive-(negative-) drift cases. 

5. In galaxies with SBHs less massive than the Milky Way’s, core formation by resonant relaxation is likely to be 
progressively less affected by anomalous relaxation. 

A. Hamers kindly provided data from his TP I code that were used in constraining the functional forms of the 
anomalous diffusion coefficients in lj3l I also thank him, F. Antonini and E. Vasiliev for comments that improved the 
manuscript. This work was supported by the National Science Foundation under grant no. AST 1211602 and by the 
National Aeronautics and Space Administration under grant no. NNX13AG92G. 


APPENDIX 

Properties of steady-state solutions in the anomalous-relaxation regime 


We consider solutions to the time-independent Fokker-Planck equation in 7?.-space in the anomalous-relaxation 
regime. We assume that the diffusion coefficients are modified versions of the resonant diffusion coefficients: 

(A7^)=n;l(£,7^)(A7^)RR, (Ala) 

(A7e)RR = 2A(f)(l-27e)(7i(£:,7e), {{Mlf)^Yi = ^A{E)n{l-n)g2{£,n) (Alb) 

and consider the two functional forms for {?iii,W 2 } that were considered in (JSJ a power-law modification, and an 
exponential modification. Diffusion in energy is ignored. 


(1) Power-law 

Identify wi and W 2 with gi and g 2 given respectively by Equations (|33ll and (j25p : 
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The flux coefficients, equations (ESI) , are 
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It is easy to verify that in this case, 
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In the limit 7Z <C {72.1,72.2}, these expressions become 
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The 7?.-directed flux is 


which, in the small-??, limit, becomes 
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Steady-state solutions can be characterized by either a constant, or a zero, flux. Setting tp-jz = 0 yields 


n\ 


-3A 


f(TZ) = f{TZ2){—J , 7?«{7?i,7?2}. 


(A 6 ) 
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When 7?i <C 7 ? 2 i a ^ 0 and the solution decreases toward 7? — O 5 the reverse is true when > 7 ? 2 - Setting 7?i — 7?2 
yields A = 0 and /(??) = const., the “zero-drift” solution. 

A steady-state solution with constant but nonzero flux has the small-?? form 
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for A ^ 2/3, with Cx an integration constant. For A = 2/3, 
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We compute C\ by requiring / to fall to zero at 7? = ??o (“empty loss cone”). The results are 
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'log (7?2/7?o) 

for A = 2/3. 

At most energies, 7?o <C 7 ? 2 - Assuming this inequality, / and p-jz have the following forms, depending on the value 
of A: 

1. A > 2/3, i.e. 7 ?i/7?2 > a/S. Defining p = 3A — 2 > 0, 
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The “zero-drift” case has 7?i = 7 ? 2 , A = 0, g = 2. 
3. A = 2/3, i.e. 7 ?i/ 7?2 = a/3 ' 
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It is interesting to compare the expressions for the flux to those that would obtain in the absence of anomalous 
relaxation. Repeating the analysis with gi = g 2 = 1, we find: 
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Dtz = 0, Dtztz = 2A{£)TZ{1-TZ), Ptz =-2A{£)TZ {1 - TZ) 
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R R 

Fig. 10.— Steady-state solutions in the anomalous-relaxation regime, und er two assumptions about how the resonant diffusion co¬ 
efficients are modified. Left panel: Power-law modification, equations with n = 8, 7^2 = 0.05 (shown by the vertical dashed 

line), and boundary condition /(7^) = 0 at 7^ = 0.005 (show by the vertical dotted line). 7^i/7^2 = {0^ 0.9,1,1.1,1.25}. Solid 
lines: 7^i/7^2 > 1; dot-dashed lines: 7^i/7^2 < !• Right panel: Exponential modification, equations IIA23L with 7^4 = 0.05 and 
7^3/7^4 = {0.99,0.999,0.9999,1.0001,1.001,1.01}. Solid lines: 7^3/7^4 > 1; dot-dashed lines: 7^3/7^4 < 1. 



^i/^2 ^3/^4 

Fig. 11.— Reduction in the steady-state flux due to anomalous diffusi on. Left panel: power-law modification; right panel: exponential 
modification. Both curves assume f{TV) = 0 at 7^ = 0.005, as in Figure fTOl 


SO that the steady state is characterized by 
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Again supposing that TI-q <C 7?.2, these expressions become: 
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Define a “reduction factor,” rj, as the ratio of this flux to the flux that would obtain in the presence of anomalous 
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relaxation, assuming the same value for /(7?.2). Then 
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if 7^l/7^2 > 73 
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Figures ITUl and ITT] plot more accurate expressions for f{TZ) and rj, computed using equations (IA3|) and (IA6I) . without 
assuming the smallness of TZ or TZq. 

The “zero-drift” solution has TZi = TZ 2 , A = 0 and q = 2. For these values, 


rj = 2 



(A21) 


For 7^o/7^2 = {0.5, 0.1, 0.01, 0.001}, r] « {0.35, 0.046, 9.2 x IQ-^, 1.4 x 10-^}. 

The value of TZi/TZ 2 favored in the numerical experiments was ~ 0.8, implying A ~ —0.56 and q ~ 3.7. For these 
values, 

” (i;) {^) ■ 

For 7^o/7^2 = {0.5, 0.1, 0.01, 0.001}, ry Ri {0.20,1.7 x IQ-^, 6.8 x 10"^ 2.0 x 10-i°}. 


(2) Exponential 

Next we identify {ici, W 2 } with {hi, ^ 2 } given by equations (l42ll and (|4T1) : 
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In the limit 72. <C {72.3,72.4}, these expressions become 
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In the same limit, the 72-directed flux is 
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where 6 = 
TZs > TZi, 


1 — TZ^/TZl < 0, a = 1 — TZ^/TZl > 0, hence / drops very rapidly to zero below TZ 


/(7^) ^ 


TZ 4 . Whereas for 
(A30) 


a rapid rise toward TZ — 0. This is qualitatively the same behavior as in the power-law case when (jj-ji = 0. 

In the case of a constant but nonzero flux, only the “zero-drift” solution (with TZs = 7 ^ 4 ) can be expressed in terms 
of simple functions: 
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In these expressions, Ei is the exponential function. The reduction factor defined above becomes, in this case, 
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